Superfluid turbulence from quantum Kelvin wave to classical Kolmogorov cascades 
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A novel unitary quantum lattice gas algorithm is used to simulate quantum turbulence of a BEC 
described by the Gross-Pitaevskii equation on grids up to 5760^. For the first time, an accurate 
power law scaling for the quantum Kelvin wave cascade is determined: The incompressible 

kinetic energy spectrum exhibits very distinct power law spectra in 3 ranges of fc-space: a classical 
Kolmogorov k~t spectrum at scales much greater than the individual quantum vortex cores, and a 
quantum Kelvin wave cascade spectrum on scales of order the vortex cores. In the semiclassical 
regime between these two spectra there is a pronounced steeper spectral decay, with non-universal 
exponent. The Kelvin k spectrum is very robust, even on small grids, while the Kolmogorov k 3 
spectrum becomes more and more apparent as the grids increase from 2048'^ grids to 5760'^. 
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Introduction. — Superfluid turbulence is an intriguing 
subject both in its own right for studying quantum turbu- 
lence in liquid Helium, Bose-Einstein condensates (BEC), 
and neutron stars, for example, but also for comparing to 
classical fluid turbulence, which itself is one of the grand 
challenge problems of the millenium [TJ [5]. There is a 
broadly acknowledged need for high resolution quantum 
turbulence simulations-in this Letter we strive to meet 
this with a unitary simulation of the Gross-Pitaevskii 
(GP) equation on large spatial grids up to 5760'^. 

Fundamental to superfluid turbulence is the quantum 
vortex: a topological singularity with the superfluid den- 
sity exactly zero at the vortex core 3 . fn the simplest 
case, all the quantum vortices are discrete, have the same 
strength (i.e. same quantized circulation in multiples of 
27r), and the flow is inviscid. This stands in sharp con- 
strast to classical incompressible fluid turbulence where 
the concept of a vortex is imprecise, the vortices are not 
discrete with arbitrary size and strength, and where vis- 
cosity plays an essential role. Moreover, in classical tur- 
bulence there are two strongly competing effects: sweep- 
ing of small scale vortices by large scale vortices, and 
straining of vortices by vortices of similar scales. Building 
on the idea of Richardson's local cascade of energy from 
large to smaller and smaller vortices till viscosity dissi- 
pates the smallest ones into heat [4], Kolmogorov [5] as- 
sumed there is an inertial energy spectrum that depends 
only on the energy input and wave number. Assuming 
the energy transfer and the interacting scales are purely 
local and sweeping is not important, Kolmogorov derived 
the inertial energy spectrum for classical incompressible 
turbulence: E{k) — Ck£^ k^i , for some constant Ck- 

Phenomenology of quantum turbulence. — Quantum 
turbulence is envisaged to arise from dense quantum vor- 
tex tangles [6] and this is borne out by numerical sim- 
ulation, for example as shown in Fig. [T] Since the flow 



outside a quantum vortex core is simple potential flow, 
it is thought that for scales ^> ^ (the radius of a quan- 
tum vortex core) the discrete nature of the quantum vor- 
tices is lost and the superfluid density is approximately 
constant. Since sweeping is now much less of an issue, 
large scale quantum turbulence could resemble incom- 
pressible classical turbulence with Kolmogorov energy 
cascade E{k) « fc~3, for k <C C^^. 

For length scales on the order of the vortex core one 
needs to consider the effects of vortex reconnection — a 
reconnection that occurs in superfluids without the need 
of any dissipative mechanism, unlike classical vortex re- 
connection processes. During the quantum vortex-vortex 
reconnection/coUision, the vortex lines are sharply dis- 
torted and emit large amplitude Kelvin waves (large rel- 
ative to the wavelength) . The Kelvin waves now interact 
with each other to generate Kelvin waves of smaller and 
still smaller wavelength. This Kelvin wave energy cas- 
cade continues until one reaches such short scales that 
phonons are emitted [TJ |2] . Just as the dissipation wave 
number fcdiss = £t- v~i cuts off the Kolmogorov ernergy 
cascade in classical turbulence, the Kelvin wave cascade 
in quantum turbulence is truncated at the wave numbers 
where sound emission breaks down. Thus a power-law 
spectrum is anticipated in the Kelvin wave cascade wave 
number range in the incompressible kinetic energy spec- 
trum E{k) Kik^°', for fc 3> ^, where the exponent a is 
yet to be determined. There has been considerable effort 
to devise theories and related numerical methods [1111117] 
that yield the precise value of this exponent as well as the 
behavior of the incompressible kinetic energy spectrum 
in the transition region between the Kolmogorov and the 
Kelvin wave cascade spectra. 

Gross-Pitaevskii equation. — At sufficiently low tem- 
peratures, the ground state wave function, Lp, of a BEC 
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FIG. 1: Quantum turbulence (zoom-in online to view tangles). 
These are vortex core isosurfaces from our quantum lattice simula- 
tion on 1024^ grid at t = 50K time step, with parameters a = 0.03 
and g = 1.004 as defined in ( |15[ |. The initial conditions are the 
same as for the 5760'^-run for a repulsive nonlinear interaction. 

can be described by the (normalized) GP equation 

idtip^ -V^ip + a{g\ip\^ ^l)ip. (1) 

We have introduced two parameters in ([T|) that are useful 
in our numerical simulations: a is simply a spatial rescal- 
ing parameter to enhance the grid resolution of the vortex 
core and 5 is a measure of the strength of the nonlinear 
coupling term in GP. Unlike the Navier-Stokes equation 
for classical turbulence, the GP equation is a Hamilto- 
nian system: the total energy Stot = const. It is well 
known [8 that the Madelung transformation ip — y/p e'^ 
on the GP equation results in compressible inviscid fluid 
equations for the density p = \(p\'^ and velocity v — 2W6, 
with the appearance of quantum pressure terms in the 
momentum and energy equations. The -Etot can be split 
into an incompressible and compressible kinetic energies, 
an internal energy and a quantum energy: 

-Btot = El°nt) + E~{t) + E^t) + = const. 

(2) 

Typically, the GP equation has been solved numeri- 
cally by the split Fourier time method 9 . In the seminal 
work of Nore et. al. [TU] , the GP equation was solved on a 
512^ grid and their incompressible kinetic energy spectra 
while not incompatible with the Kolmogorov k~'s scaling 
did not unequivocally prove the k~ 3 scaling. Barenghi [2J 
and Kobayashi et. al [1] attributed this to the presence 
of the Kelvin wave cascade. The Tsubota group [1] in- 
troduced a wave number dependent dissipative term into 
their simulations to damp out wave numbers on the or- 
der of the vortex core. While this suppresses the Kelvin 
wave cascade on the quantum turbulence, it also leads 
to a time decay in both the total number and total en- 
ergy i?TOT- To circumvent the decay in the total number, 
they add a time-varying chemical potential in the GP 



equation although the total energy still decays. Most 
of their simulations were restricted to a 512'^ grid and 
did not yield a convincing incompressible kinetic energy 
spectrum of k~3 for this augmented GP equation. These 
earlier simulations could not (nor did not want) to resolve 
the Kelvin wave cascade regime and its spectral power. 

In this Letter we will, for the first time (to our knowl- 
edge), perform very high grid resolution runs (up to 
5760'^) of the Hamiltonian GP equation to resolve the 
quantum Kelvin wave cascade. We use a novel unitary 
quantum algorithm whose solution clearly exhibits three 
power law regions for E'j^^°°"'^{k): for small k the Kol- 
mogorov fc~3 spectrum while for high k a Kevlin wave 
spectrum of k~^. A transitional power law on the order 
of k^^ to fc^^ joins these two spectral regions. For very 
large k, the phonon emissions abruptly cut off the Kelvin 
wave power law spectrum. 

Unitary quantum lattice gas algorithm. — To recover 
the GP equation in 3-1-1 dimensions we will need only 
specify 2 complex numbers per point a; on a cubic lattice 
and consider 

^(-'^)=(?&1)) 

where in a quantum computer a and f3 would be the "ex- 
cited" state amplitudes of two qubits , respectively, at x. 
At each point, these excited state probability amplitudes 
are entangled by a collision operator of the form 

gif<T.(i-^.)^ (4) 
where the Pauli spin matrices are 

^ (1 0) = (z 0*) = (0 -1) ■ 

(5) 

With n = ^(1 — (Jz) and n — ^(1 + cTz), the local qubit 
entanglement is then spread throughout the lattice by 
streaming operators 

^Ax,o = n + e^"^- n, 5ax,i + e^"^- n, (6) 

unitarily shifting the components of ip along ±Ax, re- 
spectively. In particular, let us consider the evolution 
operator for the 7th component of -0. Our quantum al- 
gorithm interleave the noncommuting collide and stream 
operators, i.e. [S^x.-yjC] ^ 0, 

^x^ = S^j^x.-yC^ S^x,^^ (7) 

where is the adjoint of C and 7 is either or 1 corre- 
sponding to streaming either the a or /3 component of ijj 
in ([3|. Since |Aa;| is small, (|7| is close to unity. 

Consider the following evolution operator for the 7 
component of ip 

C/,[f}(=r)] = 4%/2^42^e-^"(-), (8) 
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where e is a small perturbation parameter and ft will be 
specified later. With this evolution operator, the time 
advancement of the state ip is given by 



ip{x,t + At) ^U^ [n] il){x,t). 



(9) 



After considerable algebra, it can be shown that for the 
particular unitary coUide-stream operators in Q and ^ , 
one obtains the following quantum lattice gas equation 
on expanding in e 



iIj{x, t + At) = %jj{x, t) - ie^ 

(-1)T£3 



'ipix,t) 



K + a,)V^V(a=,t) + 0(e*), 



(10) 



where 7 = or 1. Since the order term in ( 10 1 changes 
sign with 7, one can eliminate this term by using a sym- 
metrized evolution operator 



u[n] = f/i 
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(11) 



Under diffusion ordering, in the scaling limit 
[tlj{x,t + At) — ^{x,t)] — > e^dttp{x,t), the quantum 
map tpix, t + At) — U[n{x)] ^p{x, t) leads to a represen- 
tation of the two-component parabolic equation 
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ijix,t) + Oie^), (12) 



where we still have not specified the operator fl. To re- 
cover the GP equation ([T]) , one simply rescales the spatial 
grid V a~^\/, contracts the 2-component field to the 
(scalar) BEG wave function ip 

^=(l,l)-^ = a + /3 (13) 

and chooses fl = .g|<y3p — 1 : 

tdtip = -VV + a{g\^\^ - 1)^ + ©(e^). (14) 

Several comments arc in order here. 1 Our unitary 
algorithm totally respects the Hamiltonian nature of the 
GP equation. No artifical numerical dissipation is in- 
troduced by our mesoscopically reversible algorithm. 2 
While the quantum lattice algorithm is ostensibly second 
order accurate, it turns out that the actual accuracy of 
our method approaches pseudo-spectral accuracies. This 
enhanced accuracy has also been noted in a somewhat 
related mesoscopic algorithmic cousin — the lattice Boltz- 
mann scheme [TTIIT^IT^. It can be traced to the imposi- 
tion of detailed-balanced collisions in the mesoscopic col- 
lision routine and emergent Onsager reciprocity. 3 This 
basic quantum lattice gas code with interleaved unitary 
coUide-stream sequence and operator fl has been bench- 
marked against exact collisions of scalar and of vector 



soliton solutions of the ID and 2D Nonlinear Schrodinger 
equation [Ml [15] as well as that for the Korteweg de Vries 
(KdV) solitons [16]. 4 Since our quantum algorithm con- 
sists of purely local collisions and streaming of informa- 
tion to the nearby grid points, it is ideally parallel. In 
fact, we have seen no saturation in performance up to 
the maximum number of cores available to us: 12,288 
cores on the GRAY XT-5 and 131,072 cores on the IGM 
Blue Gene/P. 5 While not stressed here, our quantum 
lattice-gas algorithms can also be run on quantum com- 
puters when they become available; two qubits are used 
to represent the -0 field at a point instead of just two com- 
plex amplitudes in ^ and the collision operator ^ is 
represented by a 2-qubit -^/swAP quantum-logic gate that 
creates pairwise entanglement between the on-site qubits. 
This entanglement spreads through the qubit system by 

Quantum turbulence simulations. — Most of our simu- 
lations had as initial conditions a set of 12 straight line 
vortices consisting of three groups of 4 vortices, with the 
group axes in the x, y and z directions. Because the 
space is periodic, these lines are topologically loops. The 
groupings by 4 was to ensure periodicity. Asymptotically, 
one can determine the form of a straight line vortex with 
unit winding number using a Pade approximate to the 
steady state solution of ( 14 1 , following Berloff [T7] . In 
polar coordinates (r, 0, z), one such straight line 0- vortex 
(centered at the origin) is 



tp{r) 



' Ilar2(12 + ar2) 
.g[384 + ar2(128+ llar2)]' 



(15) 



asymptotically. \ip\ — > l/y^ as r — > 00, and \ip\ ^ aj g 
as r — > 0. A typical isolated core radius scales as 
^ ~ (8/a) 2 . Two or more perpendicularly oriented line 
vortices are unstable and vortex entanglement ensues 
from such initial conditions. 

The incompressible kinetic energy spectrum can be ex- 
tracted from the conserved total energy of the system, 
following Nore et. al [5]. On the top of Fig. |2]we present 
such spectra from our simulation of the GP equation on 
a 2048"^ grid for a = 0.02 and g = 3 at evolution time 
t = 8i^ and t = 20if (in lattice units). The power laws 
are determined by linear regression. Thus, within a sin- 
gle simulation run, we find that the incompressible ki- 
netic energy spectrum has three distinct power law fc~" 
regions that range from the classical turbulent regime of 
Kolmogorov for "large" scales (much greater than the in- 
dividual quantized vortex cores) to the quantum Kelvin 
wave cascades at the "small" scales (on the order of the 
individual quantized cores). There is a semi-classical 
region adjoining the Kolmogorov and Kelvin spectra. 
These three power-law regions are quite robust as shown 
in middle and bottom of Figs. [2] from simulations on 
larger grids: 3072'^ and 5760^ and different initial con- 
ditions. In particular, the initial conditions for the 5760^ 
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FIG. 2: The incompressible kinetic energy spectra for a periodic 
12-vortex set with a = 0.02, and an initial core radius g 15 it 5. 
The linear regression fits for power-law k~" yield a's given in Table 
I. There are 3 distinct spectral regions: (a) k~3 Kolmogorov energy 
cascade for small k, (b) steep semi-classical transition region for 
intermediate k, and (c) k~^ Kelvin wave cascade for large k. The 
Kolmogorov cascade becomes robust for large grids, as seen by the 
insets. 



grid simulation are chosen to have very long Poincare 
recurrence time. Since the GP equation is Hamiltonian, 
Poincare recurrence exists for arbitrary initial conditions. 
We found Poincare recurrence to occur very rapidly for 
these simple 12-vortex systems (because of space limita- 
tions we shall discuss these results elsewhere). Hence to 
have very large Poincare recurrence times we chose initial 
conditions of the form <i> = (11^1^9^)^. The exponent a for 



the spectral are given in Table |l] together with the 
range of k for these regions. A linear regression fit for the 
Kolmorogov range is shown in Fig. [2] A similar fit for the 
Kelvin wave range is not shown since there is no point 
scatter about the regression line. The sharp drop-off in 
the spectrum at the end of the Kelvin wave cascade is 
due to the emission from very short wavelength phonons. 



TABLE I: The spectral exponent k " for the 3 distinct regions, 
as determined by linear regression. The first row for 2048^ grid is 
at time t = 8K while the next row for 2048^grid is at time t = 20K 



Grid 


Kolmogorov 


Semi-classical 


Kelvin wave 


2048^ 


1.73 (6 < fc < 30) 


6.59 (60 < fc < 140) 


2.96 (250 < fc < 600) 


2048^ 


1.84 (6 < fc < 30) 


6.34 (60 < fc < 140) 


2.97 (250 < fc < 600) 


3072=* 


1.69 (7 < fc < 45) 


7.11 (120 < fc < 200) 


3.01 (220 < fc < 1000) 


5760=* 


1.68 ( 90 < fc < 230) 


7.12 (430 < fc < 600) 


3.00 (1000 < fc < 1650) 
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